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Current-voltage characteristics and the linear resistance of the two-dimensional XY model with 
and without external uniform current driving are studied by Monte Carlo simulations. We apply 
the standard finite-size scaling analysis to get the dynamic critical exponent z at various tempera- 
tures. From the comparison with the resistively-shunted junction dynamics, it is concluded that z 
is universal in the sense that it does not depend on details of dynamics. This comparison also leads 
to the quantification of the time in the Monte Carlo dynamic simulation. 
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Equilibrium properties of the two-dimensional (2D) 
XY model have been studied extensively and there are 
well-established consensus on them.EI Dynamic properties 
of 2D XY model, however, are still under strong debate: 
For example, some studies on the dyoaxnic critical expo- 
nent z seem to contradict each other. Wu Investigation of 
dynamic properties requires a set of dynamic equations 
or impositions of update rules and there exist various 
possible choices for the XY model. On the one hand, 
there are Langevin-type stochastic equations like the re- 
sistively shunted junction (R>jU) dynamics for XY model 
with phase representation,Btinj and the relaxational dy- 
namics (often called the time-|dfl»endent Ginzburg Lan- 
dau dynamics) for both phaseSBcl and vortex representa- 
tions!] On the other hand, it is also possible to use Mcmtc 
Carlo (MC) simulations to study dynamic behaviorsB't2l 
although the identification of a Monte Carlo update step 
with a jtipe step still lacks a complete rigorous justi- 
fication.til An understanding of the dynamic transport 
properties, like the current-voltage characteristics, in the 
above dynamic models are important also in the practi- 
cal sense, since such properties are directly measured in 
experiments on Ay-like realizations such as Josephson- 
junction arrays and high-T c superconductors. 

We here use the standard MC simulation method 
to study the current-voltage characteristics of the XY 
model with phase representation. The purpose of the 
current investigation is to find the dynamic critical ex- 
ponent z in MC dynamic simulation and compare it with 
existing results from other dynamics. It then turns out 
that it is possible to quantify the MC time scales from 
the comparison with RSJ dynamics (see Ref. [l2] for a 
recent discussion on the time quantification of MC dy- 
namics). This opens a possibility to use MC simulation, 
which is usually much more efficient than direct numeri- 
cal integrations of stochastic equations of motion, to ex- 
amine many interesting long-time dynamic properties of 
the XY models. 

The fluctuating twist boundary condition (FTBCJ has 
been introduced first in the static MC simulationjij and 
later has been extended to RSJ dynamics a In the compu- 
tations of dynamic quantities, it has been shown that the 



FTBC has some advantage over the periodic boundary 
condition: For example, the fluctuating linear resistance 
can be .calculated in a much more convenient way in the 
FTBC.au We start from the RSJ dynamic equations of 
motion for the 2D XY model on a square lattice in the 
presence of an external current in the x direction, with 
the current density J = (J, 0) under the FTBC.0 The. 
phase variable 9i on the zth site at position satisfiestll 
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where Gij is the lattice Green function, the primed sum- 
mation is over nearest-neighbor sites (k) of j, A = 
(A x ,A y ) is the fluctuating twist variable, and = 
rfe — Vj. The thermal noise ?7y is white noise satisfy- 
ing (rjij) = and (rjij(t)r) k i (0)) = 2T(S ik Sji - SuS jk )S(t), 
with the temperature T in units of Josephson coupling 
strength. The equations of motion for A are written as 
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where denotes the summation over all links in 

the x direction, and the thermal noise terms satisfy 
(VA X ) = (VA y ) = (VA X VA V ) = and (■nA x ( t ) T lA x { )) = 
(VA y {t)vA v (0)) = (2T/L 2 )S(t). The above Langevin- 
type equations of motion (Q)-(^|) can be cast into the 
Fokker-Planck equation for the probability distribution 
function P: 
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whose stationary solution (dP/dt = 0) in the form 
P = epp T consequently leads to the effective Hamil- 
tonianjla 
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where the summation runs over all the nearest- 
neighboring bonds constituting the array. The Hamil- 
tonian (|^) is not invariant under the transformation 
A x — > A x + 27r, in contrast to the original RSJ equa- 
tions. However, we note that in MC simulations the en- 
ergy difference of configurations, not the energy itself, 
determines the time evolution of the system. Thus a di- 
rect use of Eq. (||) in the MC simulations does not cause 
any problem. 

In the MC simulations, we first pick one site and try to 
rotate the phase angle at this site by an amount randomly 
chosen in [—69, 69] (we call 69 the trial angle range). This 
MC try is accepted according to the standard Metropolis 
algorithm applied to the Hamiltonian (Q). After sweep- 
ing through all the lattice sites to update all the phase 
variables, we update the fluctuating twist variables A^, 
and Aj, in a similar way: we try to rotate LA X and LA y 
within the angle range 6A. (For convenience, we use 
69 = 6 A.) This complete update of the phase variables 
and the twist variables constitute one MC step. It is this 
MC step which is identified with one time step. A static 
quantity in equilibrium does not depend on the choice 
of 59; the choice of 69 only affects the convergence rate 
of the quantity (too small and too large 59 worsen the 
convergence). This is contrary to a dynamic quantity 
which in general is expected to depend on 59. In Ref. [lj, 
from the study of the anisotropic magnetic particle in 
an external magnetic field, it has been shown that as far 
as long-time behaviors are concerned the real dynamics 
and MC dynamics become in good agreements when the 
time evolution in MC simulation becomes slow enough. 
If such a result also did carry over to our model it would 
for our MC simulation imply that the MC dynamics give 
the same current-voltage characteristics as the RSJ dy- 
namics after a suitable normalization of time, when 59 is 
sufficiently small. As will be shown below, we find that 
this is indeed the case. 

In the absence of external current [J = in Eq. (g)], 
we first calculate LA X as a function of time, typically up 
to t = 10 9 -10 10 after equilibjcation, and then compute the 
fluctuating linear rcsistanceo 

Run = ^([LA x (e)-LA x (Q)} 2 ), (6) 

where 6 is chosen to be sufficiently large and (• • •) is sub- 
stituted by the time average (see Refs. || and |l6|). The 
linear resistance is closely related with the time scale r, 
defined by the time spent between the jumps of LA X by 
2tt (see Fig. 8 in Ref. §). 
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In other words, the linear resistance is obtained from the 
equilibrium fluctuations of LA X which has a time scale 
r related with jumps by 2n. Since the computation of 
i?H n in practice is more straightforward than r, we calcu- 
late the former and r is simply obtained by Eq. (0) O As 
is well known, the whole low-temperature phase of the 
2D XY model is quasi-critical and the dynamic critical 
exponent z is defined by r ~ L z due to the infinite cor- 
relation length. Alternatively, z can be computed from 
Eq. © 
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at a given temperature T. Although there exist sev- 
eral studies with a different conclusionJlj there is nev- 
ertheless a growing consensus that the 2D XY model in 
the low-temperature phase has a z which is a monoton- 
ically decreasing function of T. Figure |l| displays (a) 
i?ii n vs L at various T and (b) i?n n vs T for various 
L, for the trial angle range 59 = ir/6. As expected, 
the dynamic critical exponent z defined by Eq. (0) be- 
comes larger as T is lowered. The dependence of i?n„ 
on 69 is also studied: we find that although the value 
of i?Hn strongly depends on 59, the exponent z appears 
to saturate as 59 is decreased, as shown in Fig. |2| where 
z(59 = 7r/18) w z(tt/&) (the bigger deviations at lower 
temperatures are due to the poor convergences in i?ij n at 
59 = 7r/18 since very long time evolution is needed due to 
very large characteristic time). Since the time scale is in- 
versely proportional to the linear resistance [see Eq. (m)], 
Fig. || implies that the time scales at 59 — ir/6 and n/18 
are simply proportional to each other in a broad range 
of temperatures and system sizes; we numerically find 
t(T,L,59 = tt/6) w ct(T,L,59 = tt/18) with c w 0.14 
independent of T and L. 

If driven by an external current in the x direction, the 
twist variable A^ decreases as time goes on since the 
smaller A x gives the lower energy [see Eq. (JH])], and the 
system develops a voltage drop across the sample with 
the electric field 
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We perform simulations for L = 4, 6, 8, and 10 at 
0.65 < T < 1.5 with 59 = tt/6 and tt/18: Figure |, as an 
example, displays the current-voltage characteristics (E 
vs J) at various temperatures for L = 8 and 59 = tt/6. 
At low currents, the current- voltage characteristics be- 
come linear (E ~ J) at any temperatures in accordance 
with Ref. [l9] for RSJ dynamics. 

In the low-temperature phase of the 2D XY model, 
the dynamic scaling in Ref. rBQL for the current-voltage 
characteristics takes the formuEj 



E/J = L- z f(JL), 
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where z and the scaling function f(x) depend on T. In 
Fig. ^, the finite-size scaling is exhibited in the form 
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(E/J)L Z = f(JL) at T = 0.80 for 89 = tt/6. In this con- 
struction z is the only one free parameter and z = 3.5(1) 
is obtained. The solid line in Fig. || is obtained from R\\ n 
computed above in the absence of external currents. 

Table [| summarizes our MC results for z with and with- 
out external currents, i.e., z from Rn n without current 
driving in Fig. |^, and z from the finite-size scaling of the 
current-voltage characteristics in the presence of exter- 
nal current (see Fig. As seen from the table, these 
two different determinations are in agreement. In Ta- 
ble | we have also included the z values obtained for the 
RSJ dynamics in Ref. |], and again find agreement. From 
this we conclude that the z values from the Monte Carlo 
dynamic simulation in the present work agree with the 
corresponding values from the RSJ dynamics. 

The above agreements suggest that some features of 
the real time dynamics for the XY model can be studied 
by MC simulations. The practical advantage is that MC 
simulations are usually much more efficient than numeri- 
cal integrations of stochastic dynamic equations like RSJ 
equations. In Fig. we carry the comparison one step 
further: we directly compare the current-voltage charac- 
teristics obtained in this paper from MC at 89 = tt/6 with 
that in Ref. |l| from RSJ, for L = 8 at T = 0.90, 0.85, 
and 0.80 (from top to bottom). It is shown that multipli- 
cations of E from MC by factors 11.6, 10.9, and 10.1 for 
T = 0.90, 0.85, and 0.80, respectively, result in perfect 
agreements within a broad range of J (the temperature- 
dependence of the conversion factor has also been found 
in Ref. |l2| ) . Since the resistance is inversely proportional 
to the time scale [see Eq. (Q)] , Fig. || leads to the result 
that the time scale of the MC simulation at 89 = tt/6 for 
a given temperature is larger than RSJ dynamic simula- 
tion by the above constant factor. 

In summary, we have performed the MC dynamic sim- 
ulation of the 2D XY model subject to the fluctuating 
twist boundary condition with and without external cur- 
rents. Through the use of the standard finite-size scaling 
method, the dynamic critical exponent z has been ob- 
tained as a function of temperature, and found to be in 
accordance with RSJ dynamic results in Ref. |[ as well 
as with the scaling prediction by Minnhagen et al. in 
Ref. Ell It has also been shown that the current-voltage 
characteristics obtained from MC dynamic simulation is 
in a good agreement with that from the RSJ dynamics- 
leading to the quantification of time in MC simulation.E3 
This implies that it is possible to use MC dynamic sim- 
ulation, which is much more efficient than numerical in- 
tegration of stochastic dynamic equations in many cases, 
to study long-time behaviors of the real dynamics of the 
XY model. For example, the time-quantified MC simu- 
lation can be useful to study glassJike dynamic behaviors 
of the XY models with disorder ,Ej as well as to investi- 
gate current-voltage characteristics of the XY model in 
very low temperatures where RSJ simulation takes too 
much time to get converged results. Ej 
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TABLE I. Dynamic critical exponent z of 2D XY model 
subject to Monte Carlo dynamics. z(Rn n ) is from the scaling 
of the linear resistance at 86 = tt/6 (Fig. ^) and z(E-J) is 
from the finite-size scaling of the current- voltage characteris- 
tics (see Fig. ^). The results from the RSJ dynamics in Ref. |^ 
are also presented for comparisons. 
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FIG. 1. The linear resistance 7?n n calculated at the trial 
angle range 89 = 7r/6 in the absence of external current (a) 
vs system size L and (b) vs temperature T. In (a), the dy- 
namic critical exponent z defined as i?ij n ~ L~ z is shown to 
increase as T is decreased (see Fig. ^|for z at various 89). 
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FIG. 2. The dynamic critical exponent z from the linear 
resistance i?n n for the trial angle range 89 = n, tt/6, and 
7r/18 (from top to bottom). As 89 is decreased, z is shown to 
saturate and z{89 = n/6) coincides well with z(S9 = n/18). 
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FIG. 3. Current-voltage characteristics (electric field E 
vs current density J) for L — 8 and the trial angle range 
89 — 7r/6 at various temperatures. As J is decreased the cur- 
rent-voltage characteristics become Ohmic (E oc J) due to 
finite-size effects. 
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FIG. 4. Finite-size scaling of the current-voltage charac- 
teristics, (E/J)L Z vs JL, at T = 0.80 and 89 = tt/6. The 
dynamic critical exponent z, the only one free parameter in 
the scaling plot, is found to z = 3.5(1) at T = 0.80, in agree- 
ment with the value in Fig. ^| obtained in the absence of ex- 
ternal driving. The horizontal full line denotes Rn n L z with 
z = 3.5 obtained in the absence of an external current, con- 
firming that the simulation results with and without external 
currents are fully consistent to each other. 
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FIG. 5. Comparison of the current-voltage characteristics 
(E-J) of the MC simulation (from Fig. |) and the RSJ dy- 
namics (in Ref. |l|) for L = 8 at T = 0.90, 0.85, and 0.80 
(from top to bottom): Multiplications by factors 11.6, 10.9, 
and 10.1 for T = 0.90, 0.85, and 0.80, respectively, to E from 
MC simulations make two curves coincide very well in a broad 
range of external current density. 
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